### Replication script for: 
### "Workfare and Attitudes toward the Unemployed: New Evidence on Policy Feedback from 1990 to 2018"
### Alexander Horn, Anthony Kevins, and Kees van Kersbergen
### Comparative Political Studies

## User guide:
# Ensure all required packages are installed (in RStudio you should be prompted to install missing packages)
# Press Ctrl-Shift-S in RStudio to run this script

## NB: Replication file created in RStudio version 2023.03.0+386, with R version 4.2.3

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Load packages 

library(tidyverse)
library(foreign)
library(dplyr)
library(broom)
library(ggrepel)
library(gridExtra)
library(ggpubr)
library(lme4)
library(lmerTest)
library(effects)
library(corrplot)
library(stargazer)
library(car)
library(MASS)
library(ggplot2)
library(texreg)
library(margins)
library(ggeffects)
library(sjPlot) 
library(table1)
library(flextable)
library(kableExtra)
library(magrittr)
library(htmltools)
library(broom.mixed)
library(stringr)

# Load data
load("Workfare Replication Data.RData")

# Set colour scheme (extended fivethirtyeight theme)
ext_538 <- c("#008FD5", "#FF2700", "#A44C9B", "#FEBC00")

## Main analysis 

agg_ts <- subset(workfare, cntry != "CA" & cntry != "US")

# Values for graphing

xrange_enabling <- quantile(workfare$centred_enabling_tally, probs = c(.05, .95))
low_enabling <- quantile(workfare$centred_enabling_tally, probs = c(.05))
high_enabling <- quantile(workfare$centred_enabling_tally, probs = c(.95))

xrange_punitive <- quantile(workfare$centred_punitive_tally, probs = c(.05, .95))
low_punitive <- quantile(workfare$centred_punitive_tally, probs = c(.05))
high_punitive <- quantile(workfare$centred_punitive_tally, probs = c(.95))

# Descriptives 

obs_table_df <- workfare %>% 
    rename(Country = country) %>%
    group_by(Country, year) %>% 
    tally() %>% 
    spread(year, n) %>%
    replace(is.na(.), 0) %>%
    mutate(Total = sum(c_across(where(is.numeric)))) %>%
    ungroup() %>%
    bind_rows(summarise_all(., list(~if(is.numeric(.)) sum(.) else "Total")))

appendix_table_1 <- flextable(obs_table_df)
appendix_table_1 <- theme_vanilla(appendix_table_1)
appendix_table_1 <- bold(appendix_table_1, i = 15)

save_as_docx(appendix_table_1, 
             path = "OA1 Table 3.docx")


labels <- list(
    variables=list(preferred_strictness = "Preferred Strictness", 
                   income_decile = "Income Bracket", 
                   conservatism = "Ideology (Higher Scores = Conservative)", 
                   enabling_tally= "Enabling Tally", 
                   punitive_tally= "Punitive Tally", 
                   full_time = "Full-time Employed", 
                   part_time = "Part-time Employed", 
                   selfemp = "Self-employed", 
                   unemployed = "Unemployed", 
                   union_member = "Union Member", 
                   education_bin = "Education (Left School After 18)", 
                   male = "Male", 
                   spouse = "Has Spouse", 
                   age = "Age", 
                   baseline = "Baseline Unemployment Generosity (Level in 1980)",
                   unemp_non_standardized = "Unemployment Rate", 
                   uegen_non_standardized = "Unemployment Generosity"),
    groups=list(""))

strata <- c(list(Total=workfare))

var_table <- table1(strata, labels)
htmltools::save_html(html = var_table, 
                     file = "OA1 Table 4.html")


# Bivariate relationship

workfare$year_2 <- str_sub(workfare$year, 3, 4) 

overall_dv1 <- workfare %>%
    group_by(cntry, year_2) %>% 
    summarise(mean_dv = mean(preferred_strictness), 
              mean_workfare_tally= mean(workfare_tally)) %>% 
    unite(country_label, c("cntry", "year_2"), sep = " '") %>% 
    ggplot(aes(x = mean_workfare_tally, y = mean_dv, label = country_label)) + 
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    geom_point(size = 1, colour = alpha(ext_538[1], .75)) +
    geom_text_repel(size = 3, segment.size = 0.25, segment.alpha = .2) +
    stat_cor(colour = ext_538[3], label.x = 0, label.y = 3, size = 2.5) +
    xlab("Workfare Balance") +
    ylab("Preferred Strictness") + 
    theme(strip.background = element_rect(fill = alpha(ext_538[1], .05))) +
    guides(colour = "none", linetype = "none") 

enabling <- workfare %>%
    group_by(cntry, wave) %>% 
    summarise(mean_dv = mean(preferred_strictness), 
              mean_enabling_tally= mean(enabling_tally)) %>% 
    ggplot(aes(x = mean_enabling_tally, y = mean_dv, label = cntry)) + 
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    geom_point(size = 1, colour = alpha(ext_538[1], .75)) +
    geom_text_repel(size = 3, segment.size = 0.25, segment.alpha = .2) +
    stat_cor(colour = ext_538[3], label.x = 0, label.y = 3, size = 2.5) +
    theme(strip.background = element_rect(fill = alpha(ext_538[1], .05))) +
    guides(colour = "none", linetype = "none") +
    facet_grid(. ~ wave)

punitive <- workfare %>%
    group_by(cntry, wave) %>% 
    summarise(mean_dv = mean(preferred_strictness), 
              mean_punitive_tally= mean(punitive_tally)) %>% 
    ggplot(aes(x = mean_punitive_tally, y = mean_dv, label = cntry)) + 
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    geom_point(size = 1, colour = alpha(ext_538[1], .75)) +
    geom_text_repel(size = 3, segment.size = 0.25, segment.alpha = .2) +
    stat_cor(colour = ext_538[3], label.x = 0, label.y = 3, size = 2.5) +
    theme(strip.background = element_rect(fill = alpha(ext_538[1], .05))) +
    guides(colour = "none", linetype = "none") +
    facet_grid(. ~ wave)

components <- workfare %>%
    group_by(cntry, wave) %>% 
    summarise(mean_enabling_tally= mean(enabling_tally),
              mean_punitive_tally= mean(punitive_tally)) %>% 
    ggplot(aes(x = mean_enabling_tally, y = mean_punitive_tally, label = cntry)) + 
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    geom_point(size = 1, colour = alpha(ext_538[1], .75)) +
    geom_text_repel(size = 3, segment.size = 0.25, segment.alpha = .2) +
    stat_cor(colour = ext_538[3], label.x = 0, label.y = 3, size = 2.5) +
    theme(strip.background = element_rect(fill = alpha(ext_538[1], .05))) +
    guides(colour = "none", linetype = "none") +
    facet_grid(. ~ wave)


# Within-country trends 

# Both dvs
custom_x_labels <- c("1990", "1999/\n2000", "2008/\n2009", "2017/\n2018")
custom_y_labels <- c("2", "3", "4", "5", "6", "7", "8")

overall_ts <- agg_ts %>%
    group_by(country, wave) %>% 
    summarise(Preferred_Strictness = mean(preferred_strictness) - 5,
              Punitive = mean(punitive_tally) / -14,
              Enabling = mean(enabling_tally) / 14,) %>% 
    gather(Dependent_Variable, dv, Preferred_Strictness, Punitive, Enabling, factor_key=TRUE) %>% 
    ggplot(aes(x = wave, y = dv, colour = Dependent_Variable, linetype = Dependent_Variable)) + 
    geom_line() +
    geom_point(aes(shape = Dependent_Variable)) +
    xlab("Survey Wave") +
    ylab("Mean Preferred Strictness") +
    labs(
        colour = "Dependent Variable:",
        linetype = "Dependent Variable:",
        fill = "Dependent Variable:",
        shape = "Dependent Variable:"
    ) +
    scale_x_continuous(labels = custom_x_labels) + 
    scale_y_continuous(labels = custom_y_labels, breaks = c(-3, -2, -1, 0, 1, 2, 3),
                       sec.axis = sec_axis(~ . * 14, name = "Tally")) + 
    scale_linetype_manual(values=c("solid", "dashed", "dashed"), 
                          labels = c("Preferred Strictness", "Punitive Tallies", "Enabling Tallies")) +
    scale_colour_manual(values=c(ext_538[1], (ext_538[2]), (ext_538[3])), 
                        labels = c("Preferred Strictness", "Punitive Tallies", "Enabling Tallies")) +
    scale_shape_manual(values=c(15, 17, 19), 
                       labels = c("Preferred Strictness", "Punitive Tallies", "Enabling Tallies")) +
    theme(legend.position = "bottom", strip.background = element_rect(fill = alpha(ext_538[1], .05)),
          panel.spacing = unit(1.4, "lines")) +
    facet_wrap(. ~ country)

ggsave(file="Figure 2.png", overall_ts, width = 7, height = 7)


# Regression results 

balance_df <- workfare %>%
    unite(country_label, c("cntry", "year_2"), sep = " '") %>% 
    group_by(country_label) %>% 
    summarise(agg_balance = mean(workfare_tally))


reg_bivariate <- workfare %>%
    unite(country_label, c("cntry", "year_2"), sep = " '") %>% 
    group_by(country_label) %>% 
    mutate(mean_workfare_tally= mean(workfare_tally)) %>% 
    nest() %>%  
    mutate(model = map(data, 
                       ~lm(preferred_strictness ~ income_decile,
                           data = .x) %>% 
                           tidy)) %>% 
    unnest(model) %>%
    filter(term == 'income_decile') %>%
    mutate(row_number(),
           significance = factor(p.value < 0.05)) 

reg_bivariate <- merge(reg_bivariate, balance_df, by  = "country_label") 

reg_bivariate <- ggplot(reg_bivariate, aes(x = estimate, y = agg_balance, 
                                           colour = significance, label = country_label)) + 
    geom_point() +
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    stat_cor(colour = ext_538[3], size = 2.5)

# With all controls 
reg_full <- workfare %>%
    unite(country_label, c("country", "year_2"), sep = " '") %>% 
    group_by(country_label) %>% 
    nest() %>%  
    mutate(model = map(data, 
                       ~lm(preferred_strictness ~ income_decile +
                               full_time + part_time + selfemp + unemployed + 
                               union_member + education_bin + male + spouse + conservatism +
                               age + I(age*age),
                           data = .x) %>% 
                           tidy)) %>% 
    unnest(model) %>%
    filter(term == 'income_decile') %>%
    mutate(row_number(),
           significance = factor(p.value < 0.05)) 

reg_full <- merge(reg_full, balance_df, by  = "country_label") 

reg_full <- ggplot(reg_full, aes(x = estimate, y = agg_balance, 
                                 colour = significance, label = country_label)) + 
    geom_point() +
    geom_smooth(method = "lm", se = FALSE, 
                colour = alpha(ext_538[3], .75), linetype = "solid", 
                size = .5) +
    stat_cor(colour = ext_538[3], size = 2.5)


# Prep special labels for panels 

dropLeadingZero <- function(l){
    str_replace(l, '0(?=.)', '')
}

## Null Models

null_agg <- lme4::lmer(preferred_strictness ~ 1 +
                           + (1|country) + (1|year) + (1|countryyear),
                       data = workfare)

summary(null_agg)


ICC.Model<-function(Model.Name) {
    tau.Null<-as.numeric(lapply(summary(Model.Name)$varcor, diag))
    sigma.Null <- as.numeric(attr(summary(Model.Name)$varcor, "sc")^2)
    ICC.Null <- tau.Null/(tau.Null+sigma.Null)
    return(ICC.Null)
}

ICC.Model(null_agg)

## Empty Models (random slopes models would fail to converge here)


empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare)


# Main Analysis

disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare)    

summary(disagg)



# Graph

p <- ggpredict(disagg_no_int, terms = "centred_enabling_tally", ci.lvl = 0.835)

disagg_panel_a <- ggplot() + 
    geom_line(data = p, aes(x = x+mean(workfare$enabling_tally), y = predicted), colour = ext_538[1]) +
    geom_ribbon(data = p, aes(x = x+mean(workfare$enabling_tally), y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    geom_freqpoly(aes(x = centred_enabling_tally+mean(workfare$enabling_tally), y=(..count../sum(..count..)*5 + 5)), 
                  colour = "black", linetype = "dotted", binwidth = 2, data = workfare) +
    labs(
        title = "Panel A: Predicted Value, by Enabling Tally",
        x = "Enabling Tally",
        y = "Predicted Value"
    ) +
    coord_cartesian(xlim = c(min(xrange_enabling+mean(workfare$enabling_tally)), max(xrange_enabling+mean(workfare$enabling_tally))), 
                    ylim = c(5, 7.2), expand = FALSE) 

p <- ggpredict(disagg_no_int, terms = "centred_punitive_tally", ci.lvl = 0.835)

disagg_panel_b <- ggplot() + 
    geom_line(data = p, aes(x = x+mean(workfare$punitive_tally), y = predicted), colour = ext_538[1]) +
    geom_ribbon(data = p, aes(x = x+mean(workfare$punitive_tally), y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    geom_freqpoly(aes(x = centred_punitive_tally+mean(workfare$punitive_tally), y=(..count../sum(..count..)*5 + 5)), 
                  colour = "black", linetype = "dotted", binwidth = 2, data = workfare) +
    labs(
        title = "Panel B: Predicted Value, by Punitive Tally",
        x = "Punitive Tally",
        y = "Predicted Value"
    ) +
    coord_cartesian(xlim = c(min(xrange_punitive+mean(workfare$punitive_tally)), max(xrange_punitive+mean(workfare$punitive_tally))), 
                    ylim = c(5, 7.2), expand = FALSE) 


p1 <- margins(disagg, variables = "centred_income_decile", 
              at = list(centred_enabling_tally= seq(min(xrange_enabling), max(xrange_enabling), length.out = 10)))
p1_plot <- summary(p1)

p2 <- margins(disagg, variables = "centred_income_decile", 
              at = list(centred_punitive_tally= seq(min(xrange_punitive), max(xrange_punitive), length.out = 10)))
p2_plot <- summary(p2)


disagg_panel_c <- ggplot(data = p1_plot, 
                         aes(x = centred_enabling_tally+mean(workfare$enabling_tally), y = AME)) + 
    geom_line(colour = ext_538[1]) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    labs(
        title = "Panel C: AME of Income, by Enabling Tally",
        x = "Enabling Tally",
        y = "Average Marginal Effect"
    ) +
    scale_y_continuous(labels = dropLeadingZero) +
    coord_cartesian(xlim = c(min(xrange_enabling+mean(workfare$enabling_tally)), max(xrange_enabling+mean(workfare$enabling_tally))), 
                    ylim = c(0, 0.14), expand = FALSE) 

disagg_panel_d <- ggplot(data = p2_plot, aes(x = centred_punitive_tally+mean(workfare$punitive_tally), y = AME)) + 
    geom_line(colour = ext_538[1]) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    labs(
        title = "Panel D: AME of Income, by Punitive Tally",
        x = "Punitive Tally",
        y = "Average Marginal Effect"
    ) +
    scale_y_continuous(labels = dropLeadingZero) +
    coord_cartesian(xlim = c(min(xrange_punitive+mean(workfare$punitive_tally)), max(xrange_punitive+mean(workfare$punitive_tally))), 
                    ylim = c(0, 0.14), expand = FALSE) 


ga <- ggplotGrob(disagg_panel_a)
gb <- ggplotGrob(disagg_panel_b)
gc <- ggplotGrob(disagg_panel_c)
gd <- ggplotGrob(disagg_panel_d)

maxWidth = grid::unit.pmax(ga$widths[2:5], gb$widths[2:5], gc$widths[2:5], gd$widths[2:5])
ga$widths[2:5] <- as.list(maxWidth)
gb$widths[2:5] <- as.list(maxWidth)
gc$widths[2:5] <- as.list(maxWidth)
gd$widths[2:5] <- as.list(maxWidth)
disagg_graph <- grid.arrange(ga, gb, gc, gd,  ncol=2)

ggsave(file="Figure 3.png", disagg_graph, width = 8, height = 8)


### Prep tables

# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats

# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling Tally (Centred)", "Punitive Tally (Centred)",
                               "Income Bracket (Centred) * Enabling Tally (Centred)", 
                               "Income Bracket (Centred) * Punitive Tally (Centred)",
                               "Enabling Tally (Centred) * Punitive Tally (Centred)", 
                               "Income Bracket (Centred) * Enabling Tally (Centred) * Punitive Tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "Appendix Table 1.html")



## Robustness Check: Random slopes at level 2 and level 3

disagg_no_int_no_cntrycontrol_Full <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                                     full_time + part_time + selfemp + unemployed + 
                                                     union_member + education_bin + male + spouse + conservatism +
                                                     age + I(age*age) +
                                                     (1|country/countryyear) + (centred_income_decile-1|country/countryyear),
                                                 control=lmerControl(optimizer="bobyqa",
                                                                     optCtrl=list(maxfun=2e5)),
                                                 data=workfare)    

summary(disagg_no_int_no_cntrycontrol_Full)

disagg_no_int_Full <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                     full_time + part_time + selfemp + unemployed + 
                                     union_member + education_bin + male + spouse + conservatism +
                                     age + I(age*age) + uegen + unemp + baseline + 
                                     (1|country/countryyear) + (centred_income_decile-1|country/countryyear),
                                 control=lmerControl(optimizer="bobyqa",
                                                     optCtrl=list(maxfun=2e5)),
                                 data=workfare)    

summary(disagg_no_int_Full)

disagg_no_cntrycontrol_Full <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                              full_time + part_time + selfemp + unemployed + 
                                              union_member + education_bin + male + spouse + conservatism +
                                              age + I(age*age) + 
                                              (1|country/countryyear) + (centred_income_decile-1|country/countryyear),
                                          control=lmerControl(optimizer="bobyqa",
                                                              optCtrl=list(maxfun=2e5)),
                                          data=workfare)    

summary(disagg_no_cntrycontrol_Full)

disagg_Full <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                              full_time + part_time + selfemp + unemployed + 
                              union_member + education_bin + male + spouse + conservatism +
                              age + I(age*age) + uegen + unemp + baseline + 
                              (1|country/countryyear) + (centred_income_decile-1|country/countryyear),
                          control=lmerControl(optimizer="bobyqa",
                                              optCtrl=list(maxfun=2e5)),
                          data=workfare)    

summary(disagg_Full)



# Graph

p <- ggpredict(disagg_no_int_Full, terms = "centred_enabling_tally", ci.lvl = 0.835)

disagg_Full_panel_a <- ggplot() + 
    geom_line(data = p, aes(x = x+mean(workfare$enabling_tally), y = predicted), colour = ext_538[1]) +
    geom_ribbon(data = p, aes(x = x+mean(workfare$enabling_tally), y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    geom_freqpoly(aes(x = centred_enabling_tally+mean(workfare$enabling_tally), y=(..count../sum(..count..)*5 + 5)), 
                  colour = "black", linetype = "dotted", binwidth = 2, data = workfare) +
    labs(
        title = "Panel A: Predicted Value, by Enabling Tally",
        x = "Enabling Tally",
        y = "Predicted Value"
    ) +
    coord_cartesian(xlim = c(min(xrange_enabling+mean(workfare$enabling_tally)), max(xrange_enabling+mean(workfare$enabling_tally))), 
                    ylim = c(5, 7.2), expand = FALSE) 

p <- ggpredict(disagg_no_int_Full, terms = "centred_punitive_tally", ci.lvl = 0.835)

disagg_Full_panel_b <- ggplot() + 
    geom_line(data = p, aes(x = x+mean(workfare$punitive_tally), y = predicted), colour = ext_538[1]) +
    geom_ribbon(data = p, aes(x = x+mean(workfare$punitive_tally), y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    geom_freqpoly(aes(x = centred_punitive_tally+mean(workfare$punitive_tally), y=(..count../sum(..count..)*5 + 5)), 
                  colour = "black", linetype = "dotted", binwidth = 2, data = workfare) +
    labs(
        title = "Panel B: Predicted Value, by Punitive Tally",
        x = "Punitive Tally",
        y = "Predicted Value"
    ) +
    coord_cartesian(xlim = c(min(xrange_punitive+mean(workfare$punitive_tally)), max(xrange_punitive+mean(workfare$punitive_tally))), 
                    ylim = c(5, 7.2), expand = FALSE) 


p1_Full <- margins(disagg_Full, variables = "centred_income_decile", 
                   at = list(centred_enabling_tally= seq(min(xrange_enabling), max(xrange_enabling), length.out = 10)))
p1_Full_plot <- summary(p1_Full)

p2_Full <- margins(disagg_Full, variables = "centred_income_decile", 
                   at = list(centred_punitive_tally= seq(min(xrange_punitive), max(xrange_punitive), length.out = 10)))
p2_Full_plot <- summary(p2_Full)


disagg_Full_panel_c <- ggplot(data = p1_Full_plot, 
                              aes(x = centred_enabling_tally+mean(workfare$enabling_tally), y = AME)) + 
    geom_line(colour = ext_538[1]) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    labs(
        title = "Panel C: AME of Income, by Enabling Tally",
        x = "Enabling Tally",
        y = "Average Marginal Effect"
    ) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous(labels = dropLeadingZero) +
    coord_cartesian(xlim = c(min(xrange_enabling+mean(workfare$enabling_tally)), max(xrange_enabling+mean(workfare$enabling_tally))), 
                    ylim = c(-0.02, 0.16), expand = FALSE) 

disagg_Full_panel_d <- ggplot(data = p2_Full_plot, aes(x = centred_punitive_tally+mean(workfare$punitive_tally), y = AME)) + 
    geom_line(colour = ext_538[1]) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, fill = ext_538[1], linetype = 0) + 
    labs(
        title = "Panel D: AME of Income, by Punitive Tally",
        x = "Punitive Tally",
        y = "Average Marginal Effect"
    ) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous(labels = dropLeadingZero) +
    coord_cartesian(xlim = c(min(xrange_punitive+mean(workfare$punitive_tally)), max(xrange_punitive+mean(workfare$punitive_tally))), 
                    ylim = c(-0.02, 0.16), expand = FALSE) 


ga <- ggplotGrob(disagg_Full_panel_a)
gb <- ggplotGrob(disagg_Full_panel_b)
gc <- ggplotGrob(disagg_Full_panel_c)
gd <- ggplotGrob(disagg_Full_panel_d)

maxWidth = grid::unit.pmax(ga$widths[2:5], gb$widths[2:5], gc$widths[2:5], gd$widths[2:5])
ga$widths[2:5] <- as.list(maxWidth)
gb$widths[2:5] <- as.list(maxWidth)
gc$widths[2:5] <- as.list(maxWidth)
gd$widths[2:5] <- as.list(maxWidth)
disagg_Full_graph <- grid.arrange(ga, gb, gc, gd,  ncol=2)

ggsave(file="OA2 Figure 1.png", disagg_Full_graph, width = 8, height = 8)



### Prep tables

# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol_Full <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol_Full) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol_Full <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol_Full))

tidy_disagg_no_int_Full <- broom.mixed::tidy(disagg_no_int_Full) %>% filter(effect == "fixed")
var_disagg_no_int_Full <- as.data.frame(VarCorr(disagg_no_int_Full))

tidy_disagg_no_cntrycontrol_Full <- broom.mixed::tidy(disagg_no_cntrycontrol_Full) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol_Full <- as.data.frame(VarCorr(disagg_no_cntrycontrol_Full))

tidy_disagg_Full <- broom.mixed::tidy(disagg_Full) %>% filter(effect == "fixed")
var_disagg_Full <- as.data.frame(VarCorr(disagg_Full))

num_lvl2 <- as.numeric(sapply(ranef(disagg_Full) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg_Full) ,nrow)[2])


# Extract random effects 

var_lvl2_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "countryyear.country")$vcov, 3)
sd_lvl2_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "countryyear.country")$sdcor, 3)
var_lvl2_int_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "countryyear.country.1")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "countryyear.country.1")$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "country")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "country")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "country.1")$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "country.1")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol_Full <- signif(subset(var_disagg_no_int_no_cntrycontrol_Full, grp == "Residual")$sdcor, 3)


var_lvl2_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "countryyear.country")$vcov, 3)
sd_lvl2_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "countryyear.country")$sdcor, 3)
var_lvl2_int_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "countryyear.country.1")$vcov, 3)
sd_lvl2_int_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "countryyear.country.1")$sdcor, 3)

var_lvl3_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "country")$vcov, 3)
sd_lvl3_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "country")$sdcor, 3)
var_lvl3_int_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "country.1")$vcov, 3)
sd_lvl3_int_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "country.1")$sdcor, 3)

variance_residual_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_Full <- signif(subset(var_disagg_no_int_Full, grp == "Residual")$sdcor, 3)


var_lvl2_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "countryyear.country")$vcov, 3)
sd_lvl2_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "countryyear.country")$sdcor, 3)
var_lvl2_int_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "countryyear.country.1")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "countryyear.country.1")$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "country")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "country")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "country.1")$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "country.1")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol_Full <- signif(subset(var_disagg_no_cntrycontrol_Full, grp == "Residual")$sdcor, 3)


var_lvl2_disagg_Full <- signif(subset(var_disagg_Full, grp == "countryyear.country")$vcov, 3)
sd_lvl2_disagg_Full <- signif(subset(var_disagg_Full, grp == "countryyear.country")$sdcor, 3)
var_lvl2_int_disagg_Full <- signif(subset(var_disagg_Full, grp == "countryyear.country.1")$vcov, 3)
sd_lvl2_int_disagg_Full <- signif(subset(var_disagg_Full, grp == "countryyear.country.1")$sdcor, 3)

var_lvl3_disagg_Full <- signif(subset(var_disagg_Full, grp == "country")$vcov, 3)
sd_lvl3_disagg_Full <- signif(subset(var_disagg_Full, grp == "country")$sdcor, 3)
var_lvl3_int_disagg_Full <- signif(subset(var_disagg_Full, grp == "country.1")$vcov, 3)
sd_lvl3_int_disagg_Full <- signif(subset(var_disagg_Full, grp == "country.1")$sdcor, 3)

variance_residual_disagg_Full <- signif(subset(var_disagg_Full, grp == "Residual")$vcov, 3)
sd_residual_disagg_Full <- signif(subset(var_disagg_Full, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol_Full, 
        ~tidy_disagg_no_int_Full, ~tidy_disagg_no_cntrycontrol_Full, ~tidy_disagg_Full,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol_Full, 
        variance_residual_disagg_no_int_Full, variance_residual_disagg_no_cntrycontrol_Full, variance_residual_disagg_Full,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol_Full, 
        sd_residual_disagg_no_int_Full, sd_residual_disagg_no_cntrycontrol_Full, sd_residual_disagg_Full,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol_Full, 
        var_lvl2_int_disagg_no_int_Full, var_lvl2_int_disagg_no_cntrycontrol_Full, var_lvl2_int_disagg_Full,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol_Full, 
        sd_lvl2_int_disagg_no_int_Full, sd_lvl2_int_disagg_no_cntrycontrol_Full, sd_lvl2_int_disagg_Full,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl2_disagg_no_int_no_cntrycontrol_Full, 
        var_lvl2_disagg_no_int_Full, var_lvl2_disagg_no_cntrycontrol_Full, var_lvl2_disagg_Full,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl2_disagg_no_int_no_cntrycontrol_Full, 
        sd_lvl2_disagg_no_int_Full, sd_lvl2_disagg_no_cntrycontrol_Full, sd_lvl2_disagg_Full,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol_Full, 
        var_lvl3_int_disagg_no_int_Full, var_lvl3_int_disagg_no_cntrycontrol_Full, var_lvl3_int_disagg_Full,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol_Full, 
        sd_lvl3_int_disagg_no_int_Full, sd_lvl3_int_disagg_no_cntrycontrol_Full, sd_lvl3_int_disagg_Full,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol_Full, 
        var_lvl3_disagg_no_int_Full, var_lvl3_disagg_no_cntrycontrol_Full, var_lvl3_disagg_Full,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol_Full, 
        sd_lvl3_disagg_no_int_Full, sd_lvl3_disagg_no_cntrycontrol_Full, sd_lvl3_disagg_Full,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol_Full), nobs(disagg_no_int_Full), 
        nobs(disagg_no_cntrycontrol_Full), nobs(disagg_Full)) -> mod_stats

# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol_Full, disagg_no_int_Full, 
          disagg_no_cntrycontrol_Full, disagg_Full, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol_Full$estimate, 
                      tidy_disagg_no_int_Full$estimate, tidy_disagg_no_cntrycontrol_Full$estimate,
                      tidy_disagg_Full$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol_Full$std.error, 
                    tidy_disagg_no_int_Full$std.error, tidy_disagg_no_cntrycontrol_Full$std.error, 
                    tidy_disagg_Full$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling Tally (Centred)", "Punitive Tally (Centred)",
                               "Income Bracket (Centred) * Enabling Tally (Centred)", 
                               "Income Bracket (Centred) * Punitive Tally (Centred)",
                               "Enabling Tally (Centred) * Punitive Tally (Centred)", 
                               "Income Bracket (Centred) * Enabling Tally (Centred) * Punitive Tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA2 Table 1.html")


## Legislative styles analysis 

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_leg_style)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_leg_style)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + n_laws_total + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_leg_style)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_leg_style)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + n_laws_total + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_leg_style)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_leg_style)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_enabling_tally + centred_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + n_laws_total + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_leg_style)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_leg_style)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_enabling_tally * centred_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + n_laws_total + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_leg_style)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980", "Total Number of Laws Passed"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA2 Table 2.html")


## Simple weighting approach

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_weighted_enabling_tally * centred_weighted_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_weighted_enabling_tally + centred_weighted_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_weighted_enabling_tally + centred_weighted_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_weighted_enabling_tally * centred_weighted_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_weighted_enabling_tally * centred_weighted_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_weighted_enabling_tally + centred_weighted_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_weighted_enabling_tally + centred_weighted_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_weighted_enabling_tally * centred_weighted_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_weighted_enabling_tally * centred_weighted_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA3 Table 1.html")


## Two-year weighting approach

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_two_yr_weighted_enabling_tally * centred_two_yr_weighted_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_two_yr_weighted_enabling_tally + centred_two_yr_weighted_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_two_yr_weighted_enabling_tally + centred_two_yr_weighted_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_two_yr_weighted_enabling_tally * centred_two_yr_weighted_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_two_yr_weighted_enabling_tally * centred_two_yr_weighted_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_two_yr_weighted_enabling_tally + centred_two_yr_weighted_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_two_yr_weighted_enabling_tally + centred_two_yr_weighted_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_two_yr_weighted_enabling_tally * centred_two_yr_weighted_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_two_yr_weighted_enabling_tally * centred_two_yr_weighted_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA3 Table 2.html")


## One year lag

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus1_enabling_tally * centred_minus1_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_minus1)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus1_enabling_tally + centred_minus1_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus1)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus1_enabling_tally + centred_minus1_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus1)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus1_enabling_tally * centred_minus1_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus1)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus1_enabling_tally * centred_minus1_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus1)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus1_enabling_tally + centred_minus1_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus1)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus1_enabling_tally + centred_minus1_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus1)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus1_enabling_tally * centred_minus1_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus1)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus1_enabling_tally * centred_minus1_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus1)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA4 Table 1.html")



## Two year lag

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus2_enabling_tally * centred_minus2_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_minus2)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus2_enabling_tally + centred_minus2_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus2)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus2_enabling_tally + centred_minus2_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus2)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus2_enabling_tally * centred_minus2_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus2)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus2_enabling_tally * centred_minus2_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus2)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus2_enabling_tally + centred_minus2_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus2)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus2_enabling_tally + centred_minus2_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus2)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus2_enabling_tally * centred_minus2_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus2)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus2_enabling_tally * centred_minus2_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus2)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA4 Table 2.html")



## Three year lag

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus3_enabling_tally * centred_minus3_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_minus3)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus3_enabling_tally + centred_minus3_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus3)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus3_enabling_tally + centred_minus3_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus3)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus3_enabling_tally * centred_minus3_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus3)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus3_enabling_tally * centred_minus3_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus3)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus3_enabling_tally + centred_minus3_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_minus3)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_minus3_enabling_tally + centred_minus3_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_minus3)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus3_enabling_tally * centred_minus3_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_minus3)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_minus3_enabling_tally * centred_minus3_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_minus3)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA4 Table 3.html")



## One year lead

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls1_enabling_tally * centred_pls1_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_pls1)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls1_enabling_tally + centred_pls1_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls1)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls1_enabling_tally + centred_pls1_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls1)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls1_enabling_tally * centred_pls1_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls1)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls1_enabling_tally * centred_pls1_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls1)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls1_enabling_tally + centred_pls1_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls1)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls1_enabling_tally + centred_pls1_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls1)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls1_enabling_tally * centred_pls1_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls1)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls1_enabling_tally * centred_pls1_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls1)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA5 Table 1.html")



## Two year lead

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls2_enabling_tally * centred_pls2_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_pls2)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls2_enabling_tally + centred_pls2_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls2)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls2_enabling_tally + centred_pls2_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls2)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls2_enabling_tally * centred_pls2_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls2)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls2_enabling_tally * centred_pls2_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls2)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls2_enabling_tally + centred_pls2_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls2)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls2_enabling_tally + centred_pls2_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls2)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls2_enabling_tally * centred_pls2_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls2)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls2_enabling_tally * centred_pls2_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls2)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA5 Table 2.html")



## Three year lead

# Regressions

empty_disaggregated <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls3_enabling_tally * centred_pls3_punitive_tally +
                                      (1|country) + (1|country:countryyear),
                                  data = workfare_pls3)



disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls3_enabling_tally + centred_pls3_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls3)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls3_enabling_tally + centred_pls3_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls3)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls3_enabling_tally * centred_pls3_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls3)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls3_enabling_tally * centred_pls3_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls3)    

summary(disagg)


# Parse relevant regression data  
tidy_empty_disaggregated <- broom.mixed::tidy(empty_disaggregated) %>% filter(effect == "fixed")

tidy_disagg_no_int_no_cntrycontrol <- broom.mixed::tidy(disagg_no_int_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_int_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_int_no_cntrycontrol))

tidy_disagg_no_int <- broom.mixed::tidy(disagg_no_int) %>% filter(effect == "fixed")
var_disagg_no_int <- as.data.frame(VarCorr(disagg_no_int))

tidy_disagg_no_cntrycontrol <- broom.mixed::tidy(disagg_no_cntrycontrol) %>% filter(effect == "fixed")
var_disagg_no_cntrycontrol <- as.data.frame(VarCorr(disagg_no_cntrycontrol))

tidy_disagg <- broom.mixed::tidy(disagg) %>% filter(effect == "fixed")
var_disagg <- as.data.frame(VarCorr(disagg))

num_lvl2 <- as.numeric(sapply(ranef(disagg) ,nrow)[1])
num_lvl3 <- as.numeric(sapply(ranef(disagg) ,nrow)[2])


# Extract random effects 

var_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int_no_cntrycontrol <- signif(subset(var_disagg_no_int_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_int <- signif(subset(var_disagg_no_int, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$vcov, 3)
sd_residual_disagg_no_cntrycontrol <- signif(subset(var_disagg_no_cntrycontrol, grp == "Residual")$sdcor, 3)


var_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$vcov, 3)
sd_lvl2_int_disagg <- signif(subset(var_disagg, grp == "country:countryyear")$sdcor, 3)
var_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$vcov, 3)
sd_lvl3_int_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "(Intercept)" & is.na(var2))$sdcor, 3)

var_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$vcov, 3)
sd_lvl3_disagg <- signif(subset(var_disagg, grp == "country" & var1 == "centred_income_decile")$sdcor, 3)

variance_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$vcov, 3)
sd_residual_disagg <- signif(subset(var_disagg, grp == "Residual")$sdcor, 3)


# Combine data 
tribble(~stat, ~tidy_empty_disaggregated, ~tidy_disagg_no_int_no_cntrycontrol, 
        ~tidy_disagg_no_int, ~tidy_disagg_no_cntrycontrol, ~tidy_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Residual", NA, NA, NA, NA, NA, 
        
        "<strong>&emsp;</strong> Variance", NA, variance_residual_disagg_no_int_no_cntrycontrol, 
        variance_residual_disagg_no_int, variance_residual_disagg_no_cntrycontrol, variance_residual_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_residual_disagg_no_int_no_cntrycontrol, 
        sd_residual_disagg_no_int, sd_residual_disagg_no_cntrycontrol, sd_residual_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country-Year Clusters", num_lvl2, num_lvl2, num_lvl2, num_lvl2, num_lvl2, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl2_int_disagg_no_int_no_cntrycontrol, 
        var_lvl2_int_disagg_no_int, var_lvl2_int_disagg_no_cntrycontrol, var_lvl2_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl2_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl2_int_disagg_no_int, sd_lvl2_int_disagg_no_cntrycontrol, sd_lvl2_int_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "Country Clusters", num_lvl3, num_lvl3, num_lvl3, num_lvl3, num_lvl3, 
        
        "<strong>&emsp;</strong> Variance", NA, var_lvl3_int_disagg_no_int_no_cntrycontrol, 
        var_lvl3_int_disagg_no_int, var_lvl3_int_disagg_no_cntrycontrol, var_lvl3_int_disagg,
        "<strong>&emsp;</strong> Standard Deviation", NA, sd_lvl3_int_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_int_disagg_no_int, sd_lvl3_int_disagg_no_cntrycontrol, sd_lvl3_int_disagg,
        
        "<strong>&emsp;</strong> Variance (Income)", NA, var_lvl3_disagg_no_int_no_cntrycontrol, 
        var_lvl3_disagg_no_int, var_lvl3_disagg_no_cntrycontrol, var_lvl3_disagg,
        "<strong>&emsp;</strong> Standard Deviation (Income)", NA, sd_lvl3_disagg_no_int_no_cntrycontrol, 
        sd_lvl3_disagg_no_int, sd_lvl3_disagg_no_cntrycontrol, sd_lvl3_disagg,
        
        "", NA, NA, NA, NA, NA, 
        "N", nobs(empty_disaggregated), nobs(disagg_no_int_no_cntrycontrol), nobs(disagg_no_int), 
        nobs(disagg_no_cntrycontrol), nobs(disagg)) -> mod_stats


disagg_no_int_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls3_enabling_tally + centred_pls3_punitive_tally +
                                                full_time + part_time + selfemp + unemployed + 
                                                union_member + education_bin + male + spouse + conservatism +
                                                age + I(age*age) +
                                                (1+centred_income_decile|country) + (1|country:countryyear),
                                            control=lmerControl(optimizer="bobyqa",
                                                                optCtrl=list(maxfun=2e5)),
                                            data=workfare_pls3)    

summary(disagg_no_int_no_cntrycontrol)

disagg_no_int <- lme4::lmer(preferred_strictness ~ centred_income_decile + centred_pls3_enabling_tally + centred_pls3_punitive_tally +
                                full_time + part_time + selfemp + unemployed + 
                                union_member + education_bin + male + spouse + conservatism +
                                age + I(age*age) + uegen + unemp + baseline + 
                                (1+centred_income_decile|country) + (1|country:countryyear),
                            control=lmerControl(optimizer="bobyqa",
                                                optCtrl=list(maxfun=2e5)),
                            data=workfare_pls3)    

summary(disagg_no_int)

disagg_no_cntrycontrol <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls3_enabling_tally * centred_pls3_punitive_tally +
                                         full_time + part_time + selfemp + unemployed + 
                                         union_member + education_bin + male + spouse + conservatism +
                                         age + I(age*age) + 
                                         (1+centred_income_decile|country) + (1|country:countryyear),
                                     control=lmerControl(optimizer="bobyqa",
                                                         optCtrl=list(maxfun=2e5)),
                                     data=workfare_pls3)    

summary(disagg_no_cntrycontrol)

disagg <- lme4::lmer(preferred_strictness ~ centred_income_decile * centred_pls3_enabling_tally * centred_pls3_punitive_tally +
                         full_time + part_time + selfemp + unemployed + 
                         union_member + education_bin + male + spouse + conservatism +
                         age + I(age*age) + uegen + unemp + baseline + 
                         (1+centred_income_decile|country) + (1|country:countryyear),
                     control=lmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                     data=workfare_pls3)    

summary(disagg)


# Create table 

stargazer(empty_disaggregated, disagg_no_int_no_cntrycontrol, disagg_no_int, 
          disagg_no_cntrycontrol, disagg, 
          coef = list(tidy_empty_disaggregated$estimate, tidy_disagg_no_int_no_cntrycontrol$estimate, 
                      tidy_disagg_no_int$estimate, tidy_disagg_no_cntrycontrol$estimate,
                      tidy_disagg$estimate),
          se = list(tidy_empty_disaggregated$std.error, tidy_disagg_no_int_no_cntrycontrol$std.error, 
                    tidy_disagg_no_int$std.error, tidy_disagg_no_cntrycontrol$std.error, 
                    tidy_disagg$std.error),
          omit.table.layout = "s",
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          dep.var.labels = " ",
          dep.var.caption = "Dependent Variable: Preferred Strictness",
          column.labels = c("Baseline Model", "Without Country Controls<br> or Interaction", "Without<br>Interaction", 
                            "Without<br>Country Controls", "Full Model"),
          covariate.labels = c("Income Bracket (Centred)", 
                               "Enabling tally (Centred)", "Punitive tally (Centred)",
                               "Income Bracket (Centred) * Enabling tally (Centred)", 
                               "Income Bracket (Centred) * Punitive tally (Centred)",
                               "Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Income Bracket (Centred) * Enabling tally (Centred) * Punitive tally (Centred)", 
                               "Full-time Employed", "Part-time Employed", "Self-employed", "Unemployed",
                               "Trade-union Member", "Left School After the Age of 18", "Male", 
                               "Has Spouse", "Conservatism", "Age", "Age * Age", "Unemployment Rate (Standardized)",
                               "Unemployment Generosity (Standardized)", "Unemployment Generosity in 1980"),        
          single.row = TRUE, no.space=TRUE, 
          star.cutoffs = c(.1, .05, .01, .001), star.char = c("+", "*", "**", "***"),
          notes = c("Cells contain restricted maximum likelihood regression coefficients, with standard errors in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), notes.append = FALSE,
          type = "html", out = "OA5 Table 3.html")
